A first-principles study of oxygen vacancy pinning of domain walls in PbTiO^ 
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We have investigated the interaction of oxygen vacancies and 180° domain walls in tetragonal 
PbTiC>3 using density-functional theory. Our calculations indicate that the vacancies do have a 
lower formation energy in the domain wall than in the bulk, thereby confirming the tendency of 
these defects to migrate to, and pin, the domain walls. The pinning energies are reported for each 
of the three possible orientations of the original Ti-O-Ti bonds, and attempts to model the results 
with simple continuum models are discussed. 



PACS: 61.72.-y, 85.50.-n, 77.84.-s 
I. INTRODUCTION 

Ferroelectric materials are of intense interest for use 
in nonvolatile memory applications, in which the electric 
polarization of an array element is used to store a bit 
of information. 1,2 However, the switchable polarization 
tends to decrease after many cycles of polarization rever- 
sal during device operation, a problem that is known as 
polarization fatigue. In the last decade, polarization fa- 
tigue in ferroelectrics has been under intensive study. Al- 
though several models have been proposed to explain this 
phenomenon, 3 there is still no consensus, and many de- 
tails of the fatigue process remain unclear. Nevertheless, 
it is generally believed that defects, especially charged 
defects, play an important role. For example, a series 
of experiments has provided some understanding of how 
such defects may pin the domain walls. 4,5 In particular, 
attention has been drawn to oxygen vacancies, which are 
often the most common and mobile defects in perovskitc 
ferroelectrics. Moreover, it has been found that fatigue 
resistance can be greatly improved by replacing the Pt 
electrodes with conductive-oxide electrodes, 6 which can 
be explained in terms of the ability of oxide electrodes to 
control the concentration of oxygen vacancies in the sam- 
ple. Furthermore, the fatigue rate in Pb(Zr,Ti)C>3 films 
was found to be very sensitive to the oxygen partial pres- 
sure above the sample, suggesting that oxygen vacancies 
strongly affect the fatigue process. 7 

Mechanisms for polarization fatigue based on pinning 
of domain walls by oxygen vacancies have been dis- 
cussed phenomenologically by several authors. 3 ' 8,9 How- 
ever, in order to put such phenomenologically theories 
on a firm atomistic basis, it is important to have detailed 
first-principles calculations that can provide information 
about the structure and energetics of the ferroelectric 
domain walls, of the oxygen vacancies, and of the inter- 
actions between the two. 

While domain walls have been studied using Landau- 
type continuum theories in earlier works, 10-12 first- 
principles calculations are essential for an accurate mi- 
croscopic description of the domain walls. 13-15 Most re- 
cently, Meyer and Vanderbilt studied 180° and 90° do- 
main walls in PbTiC>3 using first-principles methods, 15 
establishing the geometry of the domain walls at the 



atomic level and calculating the creation energy of the 
domain walls. The domain wall was found to be ex- 
tremely narrow, with a width of the order of the lattice 
constant a; the positions of the atoms change rapidly in- 
side the domain wall and converge to their bulk value 
very quickly outside. As for oxygen vacancies, recent 
calculations have provided very useful information about 
the structure of oxygen-vacancy defects in this class of 
materials. 16-18 However, to our knowledge, direct first- 
principles studies of the interactions between vacancies 
and ferroelectric domain walls have not yet appeared. 

Thus, the goal of this work is to use first-principles 
calculations to examine how an oxygen vacancy would 
interact with a ferroelectric domain wall, and thus to 
shed some light on how oxygen vacancies might affect 
the switching process and cause polarization fatigue. In- 
deed, as the atomistic structure and polarization profile 
are very different in a domain wall than in the bulk of 
a ferroelectric domain, one may expect that the oxygen 
vacancies would also behave differently in such very dif- 
ferent environments. To explore these issues, we adopt 
PbTiC>3 as a model system for this study, and calculate 
the formation energies for neutral oxygen vacancies in 
the bulk and in 180° domain walls of tetragonal PbTiC>3. 
Our calculations indicate that the vacancies do have a 
lower formation energy in the domain wall than in the 
bulk, thereby confirming the tendency of these defects to 
migrate to, and pin, the domain walls. The order of mag- 
nitude of the computed pinning energy is ~ 10 -1 eV, with 
substantial variations depending on geometrical configu- 
ration. 

The rest of the paper is organized as follows. In Sec. II, 
we describe the technical details of our computational 
method and the supercells we used to model the domain 
walls. We present the results on the vacancy pinning 
energies of each of the three possible orientations of the 
original Ti-O-Ti bonds from first-principles calculations 
in Sec. III. A simple continuum model is introduced 
and used to help understand these results in Sec. IV. In 
Sec. V we discuss briefly the anticipated role of domain- 
wall pinning in fatigue, and discuss the case of charged 
vs. neutral defects. Finally, we summarize in Sec. VI. 
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FIG. 1. Schematic illustration of the 8x^/2x^/2 supercell used in the calculations, (a) View in x-z plane. For purposes 
of illustration, the cell is shown divided into "bulk regions" (polarized alternately along ±z, not shown) and "domain-wall 
regions" (comprised of the two primitive cells adjacent to the domain- wall center, which lies on a Pb-O plane), although the 
actual behavior is slightly more gradual. Vb and Vd w refer to possible locations of vacancies in bulk and domain-wall regions 
respectively, (b) View in y-z plane (with V2x\/2 supercell illustrated by the dashed lines), showing the distance between 
vacancies and their periodic images (~ a\/2). 



II. METHODOLOGY 

Our calculations are based on standard density- 
functional theory (DFT) within the local-density ap- 
proximation (LDA). We use a planewave pseudopo- 
tential code implemented by the Vienna ab-initio 
Simulation Package (VASP). 19 ' 20 Vanderbilt ultrasoft 
pseudopotentials 21 are used, with Pb 5d and Ti 3s, 3p 
electrons included explicitly. The 29 Ry cut-off used here 
was well tested and determined to be adequate for this 
material and these pseudopotentials. 15 The structure is 
considered to be relaxed when the forces are less than 
0.05 eV/A; the change of total energy at this time is typ- 
ically less than 1 meV. 

We first carried out reference calculations for domain 
walls without vacancies, following closely the approach 
of Rcf. 15. We describe the 180° domain wall using an 
8x1x1 supercell (long direction along x) with 2x4x4 k- 
point sampling; the tetragonal axis and polarization are 
along ±z, and two yz-oriented domain walls divide the 
supercell into up and down domains of equal width. 

We also carried out reference calculations for vacancies 
in the bulk. There are three kinds of oxygen vacancies 
that can be formed, depending on whether the removed 
oxygen atom had its Ti-0 bonds along x, y, or z, which 
we denote as V(x), V(y), or V(z), respectively. These 
bulk vacancies were studied using a 2 x 2 x 2 supercell with 
the tetragonal axis (and polarization) chosen along z, and 
with 2x2x2 k-point sampling. (Thus, V(x) and V(y) are 
equivalent by symmetry in the bulk.) To obtain the va- 
cancy energy, we first calculate the total energy of the 
pure tetragonal 2x2x2 supercell as the reference energy. 
One oxygen is then removed from the supercell in such 
a way that the supercell remains net neutral, and the 
resulting structure is relaxed. We keep the volume and 
shape of the supercell fixed and allow only the ion po- 
sitions to relax. The vacancy energy is then calculated 
by comparing the total energy difference before and after 
removing the oxygen atoms from the supercell. 



In order to study the interaction of the vacancies with 
the 180° domain wall, we constructed an Sxv^xv^ su- 
percell by doubling the 8x1x1 supercell in the y-z plane 
in a c2x2 or v / 2xv / 2R(45°) arrangement, and then re- 
moving one oxygen atom (and its periodic images) from 
the interior of the domain-wall structure. This is illus- 
trated schematically in Fig. 1, referring here to vacan- 
cies labeled "Vdw" (i-c, located in the domain wall) in 
Fig. 1(a). The vacancies form sheets in the y-z plane, 
with individual vacancies separated by a distance of 
about ay/2, as indicated in Fig. 1(b). Even sparser ar- 
rangements would be desirable, but the supercell already 
contains 80 atoms, and computational limitations make 
it difficult to increase the separation further. A 1x3x3 k- 
point mesh is used for these supercells, and once again the 
ion positions are allowed to relax while keeping the super- 
cell volume and shape fixed. (A single isolated domain 
wall decorated by vacancies could not expand or contract 
in the y or z directions because of the epitaxy constraint 
to the bulk, and our experience, consistent with Ref. 15, 
indicates that the lattice would expand very little in the 
x direction if allowed to do so.) 

Finally, in order to calculate the energy difference for 
an oxygen vacancy to be inside the domain wall, rela- 
tive to being in bulk, we find that it is advantageous to 
carry out corresponding 8x \/2x y/2 supercell calculations 
in which the oxygen vacancy has been removed from a 
bulk-like region of the supercell. In this way, we reduce 
the systematic errors associated with supercell shape, k- 
point sampling, etc. Our results for the binding energies 
of oxygen vacancies in the domain walls will normally be 
based on calculations of this kind, except where unavail- 
able as noted below. 



III. RESULTS 

We first report our calculations on isolated vacancies in 
bulk ferroelectric PbTiOa in the 2x2x2 supercell. Using 
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the theoretical values of the lattice constants (a = 3.86 
A, c/a — 1.0466) as obtained in Ref. 15, our calcula- 
tions show that the oxygen vacancies of type V(z) are 
more stable than V(x) and V(y) by 0.3 cV, similar to 
what was found by Park and Chadi for somewhat dif- 
ferent conditions. 18 As indicated in Sec. II, however, we 
prefer to calculate pinning energies for oxygen vacancies 
in domain walls by using reference bulk vacancy calcu- 
lations in an 8x^2x^2 supercell that can be used with 
and without domain walls, in order to provide maximum 
cancellation of systematic errors. Such calculations will 
be the basis for the the results given in this section. How- 
ever, we will return to the use of the 2x2x2 supercell in 
Sec. IV for some calculations relevant to the modeling of 
vacancies in a perturbed environment. 

We then confirm the structure of the relaxed 180° do- 
main wall. Our results for the 180° domain wall are very 
close to those of the previous calculation. 15 The domain 
wall is located on the Pb-0 plane. The atomic displace- 
ments converge rapidly to their bulk values, with most of 
the displacements ocurring within one unit cell of the do- 
main wall center. Consequently, the polarization reverses 
sharply within approximately one unit cell of the domain 
wall, and quickly saturates to a bulk value further away. 
Thus, the domain wall is extremely narrow, only about 
two lattice constant wide, and we find that it displays a 
significant x-z shear in addition to the reduced polariza- 
tion. With this justification, we can heuristically divide 
the supercell into two different regions, a bulk region and 
a domain-wall region, as sketched in oversimplified form 
in Fig. 1(a). 

To minimize the interactions between neighboring va- 
cancies, the domain-wall supercell is doubled as described 
in Sec. II. The shape of the supercell in the y-z plane 
is indicated by the dashed lines in Fig. 1(b). The to- 
tal energy E of this domain wall structure (80-atom 
Sxv^xv^ supercell) is calculated for use as the refer- 
ence energy. 

To calculate the oxygen vacancy energy, we remove one 
oxygen from the supercell and relax the resulting struc- 
ture. Recall that there are three types of oxygen vacan- 
cies, V(x), V(j/), and V(z), according to the orientation 
of the Ti-O-Ti bond from which the oxygen atom has 
been removed. For each type, we choose one vacancy as 
close as possible to the center of the bulk region, and 
another as close as possible to the central domain-wall 
plane. We label these as 'b' (bulk) and 'dw' (domain 
wall), respectively. Thus, Vd w (a;) is an x-oriented va- 
cancy at the domain wall, etc. In all calculations, we 
keep the supercell charge-neutral. 

All unrelaxed vacancies have an M y mirror symme- 
try through the defect site. In addition, the unrelaxed 
Vdw{x) and Vb(a;) defects have extra and M x sym- 
metries respectively, because they lie precisely in, or else 
half-way between, the domain-wall planes. When relax- 
ing the structure, we observe that the two Ti ions that 
neighbor the vacancy site relax most, as expected. In the 



TABLE I. Vacancy formation and pinning energies. A_E dw 
is the energy to create a vacancy in the domain-wall region 
of the supercell. AE^ is the reference energy to create a va- 
cancy either in a supercell without domain walls ( "Ref. with- 
out DW"), or in the bulk- like region of a supercell with do- 
main walls ("Ref. with DW"), and E p = AE b - AE iw is the 
corresponding pinning energy. All energies are in eV. 
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V(z) 


10.726 


10.777 


0.051 


10.764 


0.038 


V(y) 


10.943 


10.946 


0.003 


10.929 


-0.014 


V(z) 


10.939 


11.102 


0.163 


11.098 


0.159 


Relaxed 












V(aO 


10.445 


10.567 


0.122 


10.542 


0.097 


v(w) 


10.660 


10.743 


0.083 


10.736 


0.076 


V(z) 


10.459 


10.720 


0.261 







bulk region, the relaxation of Vb(x) and Vb(z) does not 
lead to any symmetry breaking, while for Vb (y) there are 
Pb ion displacements that break the M y symmetry and 
contribute about 10 meV to the relaxation energy. In 
the domain- wall region, Vdw(y) also has similar distor- 
tions, but the energy is only lowered by about 5 meV. As 
for Vdw (x) , which starts with a relatively high unrelaxed 
symmetry (four-clement point group C2fc), the relaxation 
breaks the M y and symmetries so that the only re- 
maining symmetry is inversion (point group C\h), and 
the total energy is lowered by about 20 meV. 

The total energy of the oxygen-vacancy supercell is 
then calculated for each configuration, and a vacancy 
formation energy is computed using an appropriate ref- 
erence. The results are listed in Table I. Here AEdw — 
-EvGdw + E oxy — Edw, where -E v edw is the energy of the 
supercell with the vacancy in the domain wall, E oxy is 
the energy of a free oxygen atom, and Edw is the en- 
ergy of a vacancy-free supercell with domain walls; and 
either AE b = E v+dw + E oxy - E dw ("Ref. with DW") 
or AE h = E v + E oxy - E ("Ref. without DW"), where 
-E-v+dw is the energy of a vacancy in the bulk-like region 
of a supercell containing domain walls, E v is the energy 
of a vacancy in a domain-wall-free supercell, and Eq is 
the energy of the corresponding bulk supercell containing 
neither domain walls nor vacancy. 

As shown in Table I, all the vacancies in the domain 
wall have lower formation energies than their counter- 
parts in the bulk. This indicates that there is an attrac- 
tive interaction between the vacancy and the domain wall 
and implies a positive pinning energy E p = AE-^ — AEdw 

Where available, the values for E p obtained using the 
reference supercell with domain walls is to be preferred 
(last two columns of Table I), since one expects a more 
systematic cancellation of errors in this case. However, 
in some situations this turns out not to be possible. For 
the vacancy V(z), in particular, a problem arises. We 
find that when we attempt to compute the relaxed en- 
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ergy E v+( \w of the supercell in which the vacancy has 
been placed as far from the domain walls as possible, 
the domain wall actually shifts its position during the 
relaxation in order to coincide with the vacancy, thus 
spontaneously converting the supercell from £v+dw to 
E ve dw This will be discussed further in Sec. Ill A. For 
this case, our best value for E p of 261 meV is obtained 
by falling back to the use of a reference supercell without 
domain walls (middle columns of Table I). By looking 
at other cases (i.e., V(x), V(y), and unrelaxed cases), it 
can be seen that an uncertainty of approximately 25 meV 
is introduced by the use of the less preferable reference 
supercell. 

It may be noted in Table I that the formation ener- 
gies AEb are significantly different for Vb(x) and Vb(y) 
(^200 meV), although by symmetry they should be equal 
in a true bulk environment. Here the differences arise 
from a supercell size effect connected with the arrange- 
ment of vacancies into sheets of fairly high density in the 
y—z plane. For a sufficiently large supercell we would 
expect these energies to become equal, because the lo- 
cal environments of Vb{ x ) and Vb(y) would be almost 
identical and the interactions between them would be 
negligible. However, in our case the vacancies are only 
separated by a distance of about a-\/2, so the interactions 
arc not negligible. On the other hand, the vacancies in 
the domain walls should have similar interactions, and 
we can expect some cancellation of errors when arriving 
at the pinning energy. Thus, we have more confidence in 
the E p values than in the relative formation energies of 
V(x), V(y) and V(z). 



A. Domain-wall shift 

As indicated in the previous subsection, when we at- 
tempt to calculate the energy E v+ dw of the Vb(z) va- 
cancy in the bulk- like region of an 8XV2XV2 supercell 
containing domain walls, the domain wall spontaneously 
shifts toward the vacancy during the relaxation process. 
We carried out further tests using a lOxv^xv^ super- 
cell and observed the same phenomenon, as shown in 
Fig. 2. First, we put the vacancy on a Ti0 2 plane inside 
the domain wall (x — 4.5a) and allow relaxation. We 
observe that the Pb-centered domain wall shifts towards 
the vacancy and becomes a Ti-centered domain wall. We 
then put the vacancy on a TiC>2 plane near the domain 
wall (x = 3.5a), and observe that the domain wall center 
(originally at x = 5a) again shifts to the left and becomes 
centered roughly at the TiC>2 plane at x — 3.5a, ending 
up with almost same total energy as before. When we 
attempt to put the vacancy even farther from the domain 
wall at x — 2.5a, we find that a new domain wall forms 
at the vacancy position. Clearly, the domain wall is try- 
ing to shift its position in each case so as to minimize 
the polarization at the position of the vacancy, thereby 
demonstrating directly the pinning effect of oxygen va- 
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FIG. 2. Circles: layer-by-layer Pb-atom z coordinates for a 
relaxed 10xv2x\/2 supercell with domain walls centered at 
x = and x = 5a. Diamonds: same, but with a layer of V(«) 
vacancies inserted at x = 4.5a. Squares: same, but with the 
vacancy layer at x = 3.5a. The domain wall can be seen to 
shift toward the vacancies. 



cancies and ferroelectric domain walls. 

This effect is most pronounced for the case of Vb(z) 
because it has the strongest pinning energy, as can be 
seen in Table I. In the case of V(y), we find a similar but 
weaker effect. That is, if the V(y) vacancy is placed close 
enough to the domain wall (e.g., at x = 4.5a in Fig. 2), 
a similar shift can occur; but no shift occurs if the defect 
is placed farther away. 

To pursue the calculation of -E v +dw in order to obtain 
AEb for Vb(z), it would be necessary for us to use a 
supercell larger than ^/2x^/2 in the y-z plane. However, 
as this would be computationally prohibitive, we have 
instead chosen to recalculate the bulk vacancy energies 
using an 8x supercell without domain walls, as 

indicated in the previous subsection. This provides an 
alternative reference energy which, though less accurate, 
is available in all cases. 



IV. MODEL CALCULATION 

Our first-principles calculations give us a rough picture 
of the interactions between oxygen vacancies and domain 
walls. These calculations show that the domain walls can 
indeed be pinned by oxygen vacancies. We obtain esti- 
mates for the pinning energies, and find that the pinning 
energy for V(z) is much larger than for V(x) or V(y). 
We would like to understand better the physics underly- 
ing these results, and to appreciate which results might 
generalize to other situations (e.g., other ferroelectrics 
materials, or other domain-wall structures such as 90° 
boundaries) . 
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With this motivation, we consider a model in which 
the vacancy formation energy depends on the immediate 
vacancy environment as characterized by local polariza- 
tions and strains. In particular, we consider a continuum 
description of the polarization and strain fields in the 
domain wall, as would occur in a Landau-type model, 
and assume that the vacancy energy can be expressed 
as a function of the local strain and polarization only. 
(A more sophisticated model might involve also a depen- 
dence on the local gradients of these fields, but we have 
not pursued this here.) Thus, in general we would write 
the vacancy formation energy as E v (r),P), where 77 and 
P are the strain tensor and polarization vector describing 
the state of the local environment before removal of the 
oxygen atom. However, we specialize here to the case 
of interest, a 180° wall lying in a y-z plane separating 
tetragonal phases with polarizations along ±z on either 
side. Thus, we focus on only the z component of po- 
larization, and for convenience we define a dimensionless 
reduced polarization p = P z / Pbuik- As for the strain ten- 
sor, we have t] yy = r\ yz = rj zz — by lattice continuity 
and r\ xy — by symmetry. Moreover, we find that r\ xx 
remains quite small in the domain- wall region. Thus, 
we focus only on the xz shear strain component and let 
i] s = (a/c)i] xz . The vacancy formation energy E v (r] s ,p) 
is thus considered as a function of the local r] s and p. 

Expanding E(-q Sl p) in powers of p, 



E v (r, s ,p) = A( Vs ) + B(r, s )p 2 + 0(p 4 ) . 



(1) 



where odd powers in p have been dropped by symmetry, 
and only terms up to 0(p 2 ) are retained henceforth. We 
then expand A(rj s ) and B(r] s ) in powers of rj s as 



and 



A(r] s ) = a + afq s + a 2 i] 2 + 0{rf s 



B(r ls )=b Q + b 1 r ls +b 2 r ] 2 +0(r 1 3 s ). 



(2) 



(3) 



Dropping terms beyond quadratic order in i] S} we take 
the coefficients a\, a 2 , 03, 61, b 2 and 63, in Eqs. (1-3) to 
constitute the parameters of our model. 

To obtain these six parameters, we do a series of cal- 
culations at a set of different rj s values for both p = 
and p = 1. To do this, we construct 2x2x2 supercells at 
different fixed values of r/ s and calculate the vacancy en- 
ergies as in the last section, allowing relaxation of ionic 
positions but not strains. A 2x2x2 k-mesh is used in 
all these calculations. We do calculations at p — by 
enforcing inversion symmetry of the lattice. If there is 
no constraint of symmetry imposed, the ions relax freely, 
resulting in p = 1. The results of these calculations are 
shown in Fig. 3. 

It might be thought that the coefficients a\ and 61 in 
Eqs. (2-3) should vanish by symmetry, but if the oxy- 
gen vacancy induces a distortion which lowers the lattice 
symmetry as discussed in Sec. Ill, this may not always be 
true. Consider, for example, the case of V(x) at -q s = 
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FIG. 3. Symbols: calculated vacancy formation energies vs. 
shear strain as obtained from 2x2x2 supercell calculations for 
(a) V(a;), (b) V(j/), and (c) V(z). Plus signs and crosses are 
for p — 1 and p — respectively. Lines are fits to Eqs. (1-3). 



and p = (inversion symmetry imposed). While the un- 
relaxed defect has T) 2 h symmetry, the relaxation lowers 
the symmetry to C 2 h (E, I, M y and Cf ) and reduces the 
energy by about 65 meV relative to the case where no 
symmetry breaking is allowed. Actually, there are two 
equivalent defects, related to each other via the broken 
symmetry M x , having degenerate energies at r/ s = but 
having a\ coefficients of opposite sign (i.e., opposite re- 
sponse to an applied xz strain r] s .) In Fig. (3) only the 
energy of the more stable of these two defects is plotted 
as a function of r) s > 0, and the dashed line is a fit to 
Eq. (2). The resulting (negative) value of a\ is given in 
Table II. 

Considering the other vacancies at p = 0, vacancy V(y) 
has a similar symmetry breaking but its orientation is 
such that the degeneracy would be split by a r\ yz strain, 



•5 



TABLE II. Coefficients obtained by fitting to the model 
calculations. (For V(x), b\ and 62 are not needed for the 
pinning energy and thus are not reported.) 





do 


ai 


«2 


bo 


61 


b 2 


V(x) 


10.705 


-1.434 


-33.117 


0.032 






v(w) 


10.705 





-8.737 


0.032 





12.315 


V(z) 


10.389 





-32.008 


0.462 


-0.994 


18.307 



which is absent here. For V(z) we observe no symmetry 
breaking at p = 0. Thus a\ vanishes for these cases. 

Turning now to p = 1, we find a reversed situation: 
V(x) and V(y) show no breaking of their G 2v symmetry 
at 7] s — 0, while V(z) breaks from to C\h after re- 
laxation. Thus, 61 = for V(y) but not for V(z). For 
non-zero rj s in Fig. (3), no further symmetry breaking is 
observed beyond what was already present at rj s = 0. 

The parameters resulting from all the fits of Eqs. (1-3) 
are listed in Table II. From the fact that ai and a 2 are 
negative, we see that the vacancies will prefer an envi- 
ronment of high shear strain. Similarly, since bo is posi- 
tive, they will prefer an environment of low polarization. 
Thus, the parameters are suggestive of a tendency for the 
vacancies to pin the domain walls. 

In order to model quantitatively the vacancy formation 
energy in the domain wall, we now have to estimate the 
values of p and r/ s that occur in a vacancy-free domain 
wall at the location where the vacancy would occur. Since 
Vdw(a^) lies in the Pb-0 symmetry plane, p = there. 
From the first-principles calculations of Ref. 15, we can 
estimate that p ~ 0.8 at the neighboring TiC>2 plane 
where Vd w (j/) and Vdw(z) are located. 

The estimation of i] s is more subtle. The problem is 
that the shear strain is not well defined in the domain 
walls, since it depends strongly on which of the sublat- 
tice we follow. For example, if we define the shear strain 
to be i] s = S z /c, where 5 Z , is the displacement of the ad- 
jacent atoms of same type in the z direction, we estimate 
that r? s (Pb) = 0.068, 77, (Ti) = 0.054, J? s (0(x)) = -0.029, 
»? s (0(y)) = -0.047, and ?? s (0(z)) = -0.027 in the center 
of the domain wall. This variation reflects the reversal 
of the polarization-related displacements along z as one 
passes through the domain wall along x. We could de- 
fine a mean shear as fj s — (1/5) J2i Vs(i) ~ 0.005, which 
as expected is about half of the "geometrical offset" de- 
fined in Ref. 15 (the offset occurs over approximately 
two unit cells). An alternative choice is the root- mean- 
square shear strain fj s = [(1/5) £V r] 2 (i)] 1 ^ 2 = 0.048. We 
expect that the most reasonable choice of an effective 
r]g S should lie somewhere between these two limits. For 
Vdw(^), the pinning energy we get from first-principles 
calculation is 97 meV. Comparing with Eqs. (1-3), we 
find that 77° = 0.03 gives a reasonable agreement with 
the first-principles result for this case, and we thus adopt 
this value. The shear strain should be slightly smaller at 
the location of Vd w (y) or Vdw{z), half a lattice constant 
away from the domain-wall center, but for simplicity we 



TABLE III. Environmental model parameters appearing in 
Eq. 4, and resulting pinning energy E p of the model compared 
with the best estimate from the direct DFT calculations in 
Table I. 





Vs 


P 


^modcl ( e y) 


£ P DFT (cV) 


V(aO 


0.03 


0.0 


0.105 


0.097 


V(») 


0.03 


0.8 


0.012 


0.076 


V(z) 


0.03 


0.8 


0.187 


0.261 



retain the same value of rj s = 0.03 for all three defects. 

Using the parameters from Table II and the values of 
7] s and p discussed above, we may estimate the pinning 
energy E p = E v (0, 1) - E v (r) s ,p) via 

E p = (1 - p 2 )b - (ai + b lP 2 )f] s - (02 + b 2 p 2 )r] 2 . (4) 

The resulting pinning energies are reported, and com- 
pared with the direct first-principles calculations, in Ta- 
ble III. 

Recall that the good agreement for V(x) occurs by 
construction. The agreement for V(z) is fair and that 
for V(y) is somewhat poor, although at least we have 
the correct sign of the pinning energy even in the worst 
case of V(y). Thus, we find that the present system is 
sufficiently complex that our simple model description is 
only partially successful. 

There are several reasons why this may be. As dis- 
cussed in Sec. Ill A, the domain walls may shift if the den- 
sity of oxygen vacancies is high, and the shift increases 
the pinning energy. This may help explain the under- 
estimation of the pinning energies for V(y) and V(z). 22 
Also, the use of a 2x2x2 supercell for the calculations of 
Fig. 3, from which the values of Table II were obtained, 
means that the vacancies were much closer to the dilute 
limit than was the case for the vacancy-in-domain-wall 
calculations of Table I. This is most likely the dominant 
source of the discrepancy for the case of Vd w (y) , which is 
less sensitive to rj s . Indeed, the fact that a higher density 
of oxygen vacancies in the domain wall leads to a larger 
pinning energy is consistent with a picture in which oxy- 
gen vacancies would tend to make planar clusters in the 
domain wall, thus acting to increase the pinning energy. 9 
Future tests on larger supercells in the y-z directions 
might help clarify these issues, although these still re- 
main intractable for the time being. 

Another obvious source of the discrepancies may be 
the limitations of the model. It is unsatisfying that the 
choice of r? s is so ambiguous, and it is unclear whether two 
variables (r) s and p) should suffice to describe the local 
environment. After all, the structural distortions change 
rapidly as one passes through the domain wall, so that 
it is not clear whether a Landau-type continuum model 
should be expected to capture the details of the energet- 
ics. It might be interesting to see whether a model more 
like the effective-Hamiltonian description of Bellaiche et 
ai, 23-25 which includes compositional disorder in order 
to treat alloys, could be successfully used here. 
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Nevertheless, we believe that our model captures some 
of the essential physics of the pinning mechanism of oxy- 
gen vacancies in the 180° domain walls. It gives the cor- 
rect sign and overall order of magnitude for the pinning 
energy, and correctly reflects that the pinning is strongest 
for Vd w (z), intermediate for Vdw(x), and weakest for 
Vdw(j/)- It also helps clarify the relative roles of strain 
and polarization effects in the pinning mechanism. Wc 
thus expect that it may be of some use for understanding 
other ferroelectric materials as well. 



V. DISCUSSION 

We now briefly discuss how oxygen vacancies may af- 
fect the ferroelectric switching process. As is well known, 
switching in ferroelectrics occurs not through a homoge- 
neous concerted reversal of the polarization in the bulk, 
but through the motion of domain walls separating re- 
gions of different polarization. Thus, insofar as these 
domain walls become pinned, the switching will be sup- 
pressed. 

In a pristine ferroelectric material that has a robust 
hysteresis loop and a large remanant polarization, the 
ferroelectric domain walls are presumably only weakly 
pinned by some pre-existing defects. Our calculations 
indicate that oxygen vacancies will tend to migrate into 
these domain walls over time, since they experience a 
binding energy to the domain wall of between 100 and 250 
meV. In fact, once inside the domain wall, we would ex- 
pect the vacancy to hop into the V(z) configuration, since 
this is lower in energy than the V(x) or V(y) configura- 
tions. Thus, the effective binding energy is ^250 meV, 
the value associated with the V(z) vacancy. 

If a significant number of vacancies accumulate in the 
domain wall, they in turn can act to pin the domain wall. 
When the density of such vacancies is low, they will not 
pin the domain walls strongly, and switching will still be 
able to occur. As the areal density of vacancies increases, 
however, an increasing fraction of domain walls (or in- 
creasingly large portions of individual domain walls) may 
become immobile, resulting in the decay of the switchable 
polarization. 

Of course, there are many limitations of our theoreti- 
cal analysis, and the real experimental situation could be 
much more complicated. We have studied neutral oxy- 
gen vacancies, which correspond to vacancies of charge 
+2e neutralized by electrons residing in nearby states 
of mainly Ti 3d character. In the absence of neutraliz- 
ing electrons, the vacancies may pin the domain walls 
even more strongly. (In this case, the situation becomes 
more complicated, since we can expect that the domain 
walls may acquire a tilt in order to compensate the va- 
cancy charges. That is, if the 180° domain wall is not ex- 
actly parallel to the polarization, there is a bound charge 
AP • n that can help neutralize the vacancy charges, an 
effect which may contribute to the strength of the pinning 



effect. 5 ) Alternatively, the oxygen vacancies may tend to 
aggregate into clusters 9 or to form defect complexes of 
various kinds. Finally, the polarization switching process 
is rather a complicated dynamical process that we have 
not attempted to model in detail. 

Nevertheless, we believe that our first-principles results 
serve as a first step towards understanding the possible 
role of oxygen vacancies in the pinning of ferroelectric 
domain walls. While we have investigated only one class 
of defects that may be involved in pinning, at least our 
calculations provide a lower limit for the strength of the 
pinning effect, which may be stronger if other defects 
or defect complexes play the dominant role. Our re- 
sults may also serve as input for more complex model- 
ing and simulation. For example, it could be used to 
extend a model such as that of Ahluwalia and Cao, who 
have done simulations of domain-wall formation in a 2D 
model simulation, 26 by the inclusion of vacancies into the 
simulation. Our results might also be useful in the formu- 
lation of an effective-Hamiltonian approach 27 that could 
be used to carry out finite-temperature simulations of the 
domain-wall behavior. Such studies might help quantify 
the effective strength of the pinning effect under more 
realistic conditions. 



VI. SUMMARY 

In summary, we have used first-principles density- 
functional calculations to investigate the interaction of 
oxygen vacancies and 180° domain walls in tetragonal 
PbTi03. Our calculations indicate that the vacancies do 
have a lower formation energy in the domain wall than 
in the bulk, thereby confirming the tendency of these 
defects to migrate to, and pin, the domain walls. The 
pinning energies are calculated for each of the three pos- 
sible orientations of the original Ti O Ti bonds, and are 
found to be 97 meV, 76 meV and 261 meV for V(x), 
V(y) and V(z) respectively. We also introduce a sim- 
ple continuum model with only two parameters (p, r/ s ) 
to model the results. This simple model gives pinning 
energies that agree qualitatively with the first-principles 
calculations, and we expect that it may prove useful for 
other ferroelectric systems as well. 

This work was supported by ONR grant N0014-97-1- 
0048. 
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